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Abstract 

A square-lattice model for the formation of secondary structures in proteins, the hydrogen- 
bonding model, extended to include the effects of solvent quality, is examined in the framework of 
the Bethe approximation. 

PACS numbers: 05.40.Fb, 05.20.+q, 05.50.+a, 36.20.-r,64.60.-i 



1 



I. INTRODUCTION 



There is a long tradition of using lattice models of polymers in an attempt to capture the 
essential features of the physics of polymers in solution l|, |2|, yj. The canonical lattice model 
is the self-avoiding walk (SAW) model extended by including interactions between nearest- 
neighbour visited sites on the lattice. These interactions model the quality of the solvent 
in which the polymer lies % |5j. When only one walk is considered, the model is thought to 
describe the behaviour of polymers in dilute solution. At high temperature the polymer is 
'happy' in solution, whilst at low temperature the polymer collapses and precipitates from 
solution. These two regimes are separated by the 9-transition0, Q, [sl. In what follows, we 
refer to this model as the B model. 

In the early 1990s a model was introduced to model the effects of hydrogen bonding on 
the formation of secondary structures resulting from the folding of a protein: the Hydrogen- 
Bonding self-avoiding walk, introduced by Bascle, Garel and Orlandjsl. In this model the 
presence of hydrogen bonds, essential in real proteins, was modelled by the presence of 
interactions between parallel straight sections of the walk, as shown in figure [H They studied 
the model in the Hamiltonian Walk limit, where all the sites of the lattice are visited exactly 
once. The model was e.te.ded fi.t to an de,.,tie.Q a„d >at. to anow fo. solvent effects 
by including all the interactions present in the G-point model, but with different interaction 
strengths, depending on the configuration of the walk (again see figure [T]) 

Exploring the phase diagrams of such frustrated self-avoiding walk models is not an easy 
task, in most cases requiring good quality numerical methods. This is mainly due to the 
frustration effects intrinsic to the presence of interactions between portions of the walk 
which may be arbitrarily far apart along the walk. The standard Monte-Carlo methods 
applied to these models consist in studying finite walks on the infinite lattice and examining 
the finite-size behaviour of the walk. The current favourite Monte-Carlo methods are the 



PERM methodpj, ll3|, UJ, ll5|, ll6|, the flat-PERM methodllTI] and the parallel tempering 
method 13]. This limits the use of the method to phase transitions coming from the zero- 
density high-temperature phase (the self-avoiding walk phase) but does not permit the study 
of phase transitions between different dense phases. The use of Monte Carlo in such phases 
is extremely difficult, normally requi ring the relaxation of some constraints, as is the case 
for the fluctuating-bond method 19|, l20l]. This however tends to erase the very effects we 
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wish to study. 

Another approach which has proved useful in studying such models is the use of transfer 
matrices in which the partition and correlation functions are expressed in terms of products 
of matrices enabling numerically exact calculations on infinitely long strips of finite width 10, 



211]. This method has two limitations: in practical terms the method is restricted to two 



dimensions, and the number of available widths is limited. The latter restriction may be 
alleviated by using the CTMRG method (corner transfer matrix re-normalisation group 



method), which enables calculations for much larger lattice sizes |22l. l23l |. These methods 
enable the investigation of the phase-diagram in the entire phase-space. 

The hydrogen-bonding self-avoiding walk in a solvent was recently studied in two di- 
mensions using a combination of transfer matrix and CTMRG methods HI] . It has also 
been independently studied using a modified PERM Monte-Carlo method in two and three 
dimensions {2^. 

All numerical methods are open to possible misinterpretation or artefacts, particularly 
when applied to models which include significant frustration effects. It is important to have 
some independent confirmation that the results obtained are reasonable. In this article 
we propose to provide such an independent confirmation by performing a mean-field type 
calculation in the form of the Bethe approximation for the hydrogen-bonding self-avoiding 
walk in a solvent and compare the results with the previously obtained numerical results. 

In the next section we present the model in detail. In section IIIII we apply the Bethe 
approximation to our model and in section [IV] we present our results. We finish with 
discussion and conclusions in section |Vl 



II. MODEL 



The model studied in this article involves the self-avoiding walk on the square lattice 
with interactions between non-consecutive visited nearest-neighbour sites on the lattice. 
In the standard G-point model the nearest neighbour interactions model effective interac- 
tions mediated by the solvent, given by the difference between the monomer- monomer and 
monomer-solvent affinities. These interactions are isotropic. In the current model, the inter- 
actions are split into two sets, as shown in figure [H those which specify a particular direction, 
the hydrogen bonds, and those that do not, which we shall refer to as the solvent-mediated 
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interactions. Hydrogen bonds carry an interaction energy —Sh and the others carry an 
interaction energy —e. The thermodynamic behaviour may be investigated by introducing 
the grand-canonical partition function, from which many of the relevant thermodynamic 
quantities may be calculated. The grand-canonical partition function is given by: 

Z^Y. ^"^exp {(5 {Nje + NhSh)) , (1) 

walks 

where Ni are the number of solvent-mediated interactions, and Nh are the number of hy- 
drogen bonds. The fugacity, which controls the average length of the walk, is denoted by 

K, and N is the total length of the walk. For convenience we define a = e/eh, and without 
changing the physics of the model, we may set £h = 1; this simply sets the temperature 
scale. The partition function then becomes: 

Z^Yk"" exp (/? {Nh + Nja)) , (2) 

walks 

from which we may calculate, for example, the free energy per site, /: 

Pf^-^\nZ (3) 

where fl is the number of lattice sites. The free energy as defined here is also the grand 
potential for the model in which we concentrate on the walk (rather than the lattice) and 
view the problem as grand canonical since the number of steps in the walk varies. Clearly 
the basic unit of the calculation is the lattice site, and so we choose the convention of 
referring to / as the free energy per site, which corresponds to the standard picture in Bethe 
approximation calculations. 

The fugacity K controls the average length of the walk. The average number of steps is 
given by 

im^K?^. (4) 

The average length increases as K is increased. For fixed a, if (3 is small enough, then the 
average length diverges continuously as K approaches some critical value Kc{a,P). This 
defines the self-avoiding walk line, which extends to a (half) plane as a is varied. If, on 
the other hand, (3 is large enough then the average length jumps discontinuously at some 
value of the fugacity, K — K*{a, (3). Together these two regimes define the plane Koo{ol, (3) 
on which the walk length first diverges. This plane separates the high-temperature, zero- 
density, phase from the low-temperature, dense, phases. 



(a) 



(b) 



(c) (d) 



FIG. 1: The nearest-neighbour interactions are spht into two classes, those of type (a) where four 
bonds forming two parallel lines model the hydrogen bonds, whilst the others (b, c and d) model 
the solvent-mediated interactions. Configuration (a) induces a preferred orientation, whilst the 
other configurations do not. 

III. THE BETHE APPROXIMATION 

In this section we describe briefly the Bethe approximation. For a good discussion of the 



Bethe approximaton see Ref. [25j. The model of interest is studied on the infinite Bethe 
lattice chosen to have the correct local geometry. The lattice chosen for the square lattice is 
shown in figure [21 The Bethe lattice is a hierarchical lattice built recursively from a central 
bond by adding to each extremity k new bonds. To each dangling bond we add k more 
bonds, and so on, such that no loops are formed. Due to the hierarchical nature of the 
lattice, it is possible to build up expressions for the partition function recursively. To see 
this, it is convenient to consider the lattice as being divided into two branches, left and 
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right for the example shown in figure [21 We may introduce the partial partition functions 
W^. and for the left and right-hand branch, respectively. These partition functions are 
conditional upon the state a of the central bond. In our model there are four possible states: 

1. empty (state 0), 

2. occupied with a link of the walk (state K), 

3. occupied with a 9 interaction (state 6), and 

4. occupied with a hydrogen bond (state H). 

By symmetry, the left and right branches will have the same partial partition functions, 
and so the l,r designation will be dropped. Each branch may be sub-divided into k sub- 
branches, such that the may be expressed in terms of the partial partition functions of 
the sub-branches. This procedure may be continued until the boundary bonds are reached. 
In order to do this explicitly, it is convenient to introduce the notion of the 'generation' of a 
link, n, which is simply the distance of the link from the boundary. As a concrete example, 
consider the calculation of wj^\ the partial partition function conditional on the central 
bond being occupied by a link of the walk. We must consider all the configurations on the 
bonds of the generation (n — 1), of which there are three for the 2 dimensional square lattice 
example shown in figure O which are compatible with the occupied central bond. Clearly 
there must be a bond leaving in one of the three directions, the other two bonds may be 
empty or occupied by a solvent-mediated interaction. If the bonds on the central bond 
and at generation {n — 1) line up, then the "empty" bonds may be occupied by hydrogen 
interactions. The weight is simply the sum of the Boltzmann weights corresponding 

to all these configurations, multiplied by the weight for adding the central link. To avoid 
the divergence of the partial partition functions it is convenient to introduce normalised 
partition functions w^^^ = PV^" I26I] with g„ chosen such that: 

E^i"^ = l- (5) 
This leads to recursion relations for the (normalised) partial partition functions: 

^':'' = ^j:c.,H}fK-'\ (6) 

{7,} i=l 
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FIG. 2: The Bethe lattice representation of the two-dimensional lattice. The dotted box shows the 
central bond, exhibiting the desired square-lattice geometry. 



where {ji} is the set of state of the k hnks forming generation n — 1, A^- is the Boltzmann 



states {ji} is compatible with the central state a, and zero otherwise. 

It is known that there is no phase transition on the infinite Bethe lattice, since the number 
of boundary sites grows too rapidly. However the recursion relations may be used in the 
centre of the lattice as self-consistency equations for the two point mean-field theory for the 
corresponding square lattice. In this case, we assume we have translational invariance, and 
drop the generational superscripts. The equilibrium states are then given by solutions of 
the following set of recursion relations: 



weight of the bond added at generation n, and the factor (^^,{7^ = 1 if the choice of the 




(7) 



K 

wk = — 



Wk (3(^0 + wqY + 2{wo + wq)wh + w% 



(8) 




Wk (^0 + we + Wh) 



(10) 



(9) 



q = wl + K (3{wo + wef + {2{wo + we) + wh)wh) wk 

+2e'^^K + ^e)^| (11) 



+e'^{wQ + we + wh)w 
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The partial partition functions give the contribution to one branch of the total partition 
function, the total (normalised) partition function conditioned upon the state of the central 
bond is then given by the product of the weight for the left and right branches. Each of the 
partial partition functions includes the Boltzmann weight corresponding to the state of the 
central bond, which is thus counted twice in the full partition function. This double counting 
is corrected by dividing each term by the relevant Boltzmann weight. Summing over all the 
possible states for the central bond gives the total (normalised) partition function, z: 

. = E^. (12) 

In the usual way, the probability of finding a given bond in state cr is given by the partition 
function conditioned upon this state divided by the total partition function, i.e. 

V, - '''' 

It should be noted that the density p of the walk on the lattice is simply px- 
The free energy per site may be related to z and q through the relation 

/3/= '''^""/'"°^ (14) 



for a full derivation of this expression see 26|]. For the square lattice A; = 3, hence j3f 



In^r — Ing. When multiple solutions to the recurrence relations exist, the solution with the 
lowest free energy is the stable equilibrium solution. 

IV. RESULTS 

It is instructive to see how the calculation works for the pure 0-point model. Many of 



the results in this case have already been presented by Lise, Maritan and Pelizzola 27|] using 
a different, variational, approach to the Bethe approximation. The pure B-point model is 
defined by a = 1, in which case ujh = and equation [10] is no longer needed. It is convenient 
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to recast the relations [3— El by setting 

Wk 

Xk = , 

Wo 
We 

Xe = — . 

Wo 

This leads directly to a trivial solution xk = xe = 0, corresponding to the zero-density 
phase. The other possible solutions are given by: 

x| + (3 - e'')xl + (3 - 2e^)xe + (e^ " 1) " l) = 0' (15) 



Xk 



\ 3(1 + xe) • ^''^ 

For a solution to be physically acceptable, xk and xq must be positive. Let us first consider 
the case 1 < < 3/2. Whilst K < 1/3 all the coefficients of the cubic equation (|TH|) are 
positive and there is no physically acceptable solution (one of the solutions is negative, and 
the other two are complex conjugates). When = 1/3 there is a solution xq = and 
Xk = 0, and when K > 1/3 this solution has a negative free energy. This is shown in 
figure El This corresponds to the high-temperature transition, where the transition point 
corresponds to an infinite self-avoiding walk, defining Koo{P) for [3 < [3q. 

When = 3/2, the coefficient of the xq term vanishes, and when = 1/3 the zero- 
density solution becomes a double root of the cubic equation (|T5l) . This change of behaviour 
is identified with the G-point, and the value /5e = ln(3/2) = 0.405465 ■ ■ ■ agrees with the 



value given by Lise et a/[27|. When > 3/2, a physically acceptable solution now exists 
for values of ii" < 1/3. To find the transition line for > 3/2 we must check for the 
stability of this new solution; for small enough K the free energy is positive, and so the 
solution corresponds to a meta-stable solution, whilst for some higher value of K, f becomes 
negative. The point where / = defines the first-order low-temperature transition, defining 
KooiP) for [3 > (3q. The free energy in this case is shown in figure [Hand the density as a 
function of (3 plotted along the transition line is shown in figure [51 

The tri critical point was identified with a double root of equation ffT^ . For > 3/2 
there is a line of double roots given by: 

For (3 > (3q this gives the position where a second order transition would have occurred had 
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FIG. 5: Density, p, plotted along the Koo{(3) line for the 0-point model (a = 1). 

it not been pre-empted by the actual first order transition, and as such is identified with the 
spinodal line. 

We now extend our analysis to a < 1. The analysis follows the same lines as for the pure 
6 model; it is possible to eliminate all the parameters in terms of xe, though this is now 
the solution of a slightly more complicated equation, which may no longer be expressed in 
a simple polynomial form. The G-point extends to a line as a is varied. For (3 < (3Q{a) 
there is still a zero-density solution for this equation when K = 1/3, corresponding to 
the self-avoiding walk transition line. The G-line is again identified with a double root of 
equation (fTSll corresponding to the zero-density phase for K = 1/3, which occurs when the 
coefficient of the xq term in a small xq expansion of the equation also vanishes. This gives 
the Pe{c() through: 



For a < 1 there is the possibility of another phase: the crystalline phase, where the 
walk fills the lattice. All the bonds align with one of the lattice directions, maximising the 
number of hydrogen bonds. This phase has zero entropy (per lattice site), and its energy 
per site corresponds to the energy for one bond and one interaction. The free energy for 




(18) 
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FIG. 6: The phase diagram in the a-(3 plane, with the fugacity, K = Koo{a,f3). The soHd hne cor- 
responds to a direct first-order transition between the self-avoiding walk phase and the crystalline 
phase. The dashed line is the line of G-points, separating the self-avoiding walk phase from the 
isotropic collapsed phase, and the diamonds correspond to the first-order transition between the 
isotropic and crystalline collapsed phases. The circle shows the location of the multicritical point 
where the different transition lines meet. 

this phase is given by 

/5/cryst = -(/? + In ir). (19) 

Setting K = 1/3, and thus following the self-avoiding walk line, it is seen that /cryst = when 
/3 = ln3. If l3e{a) < In 3, the first transition met on increasing jS is the 6-point transition, 
and the location of the crystallisation transition, /^/^(a), is determined by comparing the 
free energies of the collapsed phase and /cryst along the Koo-Une. If /5e(a) > In 3, the walk 
collapses directly to the crystalline phase, and the B-transition is not present. The change- 
over between these two cases occurs at a multi-critical point, which is found by setting 
/3e(«mc) = In 3. This gives the location of the multi-critical point as: K^^ = 1/3, /^mc = 
ln3 ^ 1.0986123 and Omc = ln(21/16)/ln3 ^ 0.2475247. The phase diagram projected 
onto the Koo{ct,[3) plane is shown in figure El When [3 > Ph{c(), the i^oo-hne is given by 
K = e~^. In the high density region of the phase diagram, comparing the free energies of the 
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FIG. 7: The phase diagram in the (3 — K plane, with a = 0.5 > amc- The circle corresponds to the 
location of the ©-point transition. The transition from the p = phase to the crystalline phases is 
first order, whilst the transition from the p = phase and the Isotropic collapsed phase is second 
order for f3 < (3q. The transition between the two dense phases is first order. 

collapsed and crystalline phase in the dense region, another phase transition may be seen, 
already reported in the literature ll|, |2J] . Here this phase transition shows up as a first-order 
transition, whilst in more realistic numerical calculations there is evidence to suggest that 
it is in fact of second order ll|]. Phase diagrams for two representative cases, a = 0.5 and 
a = 0.2, are shown in figures [7] and [HI 



V. DISCUSSION 



Whilst the Bethe approximation is "only" a mean-field type calculation, it usually cap- 
tures the essential features of the model. For a = 1 we obtain results consistent with the 
results of Lise et a/R. For a = the model corresponds to the pure Hydrogej^bonding 
model. The {K,(3) phase diagram, as already remarked by Buzano and Pretti[29], corre- 
sponds closely to that found with transfer matrices lOj , except that the phase transition 
between the isotropic dense phase and the crystalline phase is found to be first order here. 
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FIG. 8: The phase diagram in the P — K plane, with a = 0.2 > amc- The transition between p = 
and the crystalhne phase and the transition between the collapsed phases are first order. The Q 
point transition is absent. 



whilst evidence suggests that it is in fact second order [10|, In the model studied here 
we include the effect of the solvent, as compared to the hydrogen type interaction, and this 
enters through the parameter a. The solvent-mediated interactions favour a collapse to an 
isotropic collapsed phase, whilst the hydrogen-bonding interaction tends to align the walk 
along one of the lattice directions, breaking the rotational symmetry, leading to a crystalline 
phase. For a close to one, the collapse of an infinite chain is progressive as j3 is increased. 
At the collapse transition the fractal dimension of the walk is less than the dimension of 
the lattice. It is expected that the details of the lattice will not influence the transition. 
The collapse transition is, in this case, in the same universality class as the standard 
point. However, once in the dense phase, the dimensions of the walk and the lattice are 
the same. The walk "sees" the lattice. In the plane where the length of walk first diverges 
{K = K*{a,(3)) we see the appearance of a second transition, from the isotropic phase to 
the crystalline phase. This transition line extends into a transition plane for K > K*{a,l3), 
as shown in figure [71 As a is lowered, a point is reached in which the hydrogen interactions 
dominate, and the walk collapses directly to a dense crystalline phase. These two regimes 
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are separated by a multicritical point, (Xmci where the three transition hnes shown in figure E] 
meet. The phase transitions found in the context of the Bethe approximation correspond 
well to the phase diagrams found numerically for the same model UJ, l2J] . Krawczyk et 
al 2J] found the transition between crystalline and isotropic dense phases to be first order 
in three dimensions, whilst in two dimensions the order of the transition was less clear, and 
the authors conjectured that the transition is critical. This conjecture is supported by the 
numerical study of Foster and Pinettes[lli]. Here this transition is found to be first order due 
to the mean-field nature Bethe approximation, which should be exact in c? = oo dimensions. 

The predicted value of amc ~ 0.248 is close to what is seen numerically (amc = 0.3 — >■ 
0.5) Uj. Similar phase diagrams to those found in figure [7| and figure [H] are found in other 
models where frustration effects are important in phase transitions between different dense 



phases, in particular the vertex-interacting self-avoiding walk|23l. |28| (figure [S]) and the bond- 
interacting walkQ, [3^ (figure [7]). However, the nature of the transition is found to be very 
sensitive to details of the interactions. The vertex interacting walk has a phase transition in 



the dense phase which is in the Ising universality class [23|, 128| whilst the Hydrogen-bonding 
model has a transition which is critical, but not Ising [y 0.87) Uj. In the context of the 
Bethe approximation, all these transitions show up as first order. It would be interesting to 
understand how to incorporate in a mean-field type calculation the essential features which 
would reproduce the second order nature of the transition between dense phases. 

The general features of what is presented here remain true in three dimensions. The 
self- avoiding walk line occurs for K = 1/5 rather than 1/3. This is easily understood: 
l/i^sAw corresponds to the average number of lattice directions available to the walk at 
each step. Due to the absence of loops on the Bethe lattice, this is simply one less than the 
co-ordination number of the lattice, i.e. 2d — 1. 
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